library(sf)
library(ggplot2)
library(leaflet)
library(scales)
library(ggmap)
library(dplyr)
ak_regions <- read_sf("../shapefiles/ak_regions_simp.shp")

plot(ak_regions) 

class(ak_regions)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
head(ak_regions)
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.2296 ymin: 51.15702 xmax: 179.8567 ymax: 71.43957
## Geodetic CRS:  WGS 84
## # A tibble: 6 x 4
##   region_id region           mgmt_area                                  geometry
##       <int> <chr>                <dbl>                        <MULTIPOLYGON [°]>
## 1         1 Aleutian Islands         3 (((-171.1345 52.44974, -171.1686 52.4174~
## 2         2 Arctic                   4 (((-139.9552 68.70597, -139.9893 68.7051~
## 3         3 Bristol Bay              3 (((-159.8745 58.62778, -159.8654 58.6137~
## 4         4 Chignik                  3 (((-155.8282 55.84638, -155.8049 55.8655~
## 5         5 Copper River             2 (((-143.8874 59.93931, -143.9165 59.9403~
## 6         6 Kodiak                   3 (((-151.9997 58.83077, -152.0358 58.8271~
ak_regions_3338 <- ak_regions %>% #use pipe because SF is part of tidyverse/dplyr
  st_transform(crs = 3338)

st_crs(ak_regions_3338)#this is a good CRS projection for AK
## Coordinate Reference System:
##   User input: EPSG:3338 
##   wkt:
## PROJCRS["NAD83 / Alaska Albers",
##     BASEGEOGCRS["NAD83",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["Alaska Albers (meters)",
##         METHOD["Albers Equal Area",
##             ID["EPSG",9822]],
##         PARAMETER["Latitude of false origin",50,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-154,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",55,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",65,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Topographic mapping (small scale)."],
##         AREA["United States (USA) - Alaska."],
##         BBOX[51.3,172.42,71.4,-129.99]],
##     ID["EPSG",3338]]
plot(ak_regions_3338)

ak_regions_3338 %>%
  select(region)
## Simple feature collection with 13 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2175328 ymin: 405653 xmax: 1579226 ymax: 2383770
## Projected CRS: NAD83 / Alaska Albers
## # A tibble: 13 x 2
##    region                                                               geometry
##    <chr>                                                      <MULTIPOLYGON [m]>
##  1 Aleutian Islands     (((-1156666 420855.1, -1159837 417990.3, -1161898 41694~
##  2 Arctic               (((571289.9 2143072, 569941.5 2142691, 569158.2 2142146~
##  3 Bristol Bay          (((-339688.6 973904.9, -339302 972297.3, -339229.2 9710~
##  4 Chignik              (((-114381.9 649966.8, -112866.8 652065.8, -108836.8 65~
##  5 Copper River         (((561012 1148301, 559393.7 1148169, 557797.7 1148492, ~
##  6 Kodiak               (((115112.5 983293, 113051.3 982825.9, 110801.3 983211.~
##  7 Kotzebue             (((-678815.3 1819519, -677555.2 1820698, -675557.8 1821~
##  8 Kuskokwim            (((-1030125 1281198, -1029858 1282333, -1028980 1284032~
##  9 Cook Inlet           (((35214.98 1002457, 36660.3 1002038, 36953.11 1001186,~
## 10 Norton Sound         (((-848357 1636692, -846510 1635203, -840513.7 1632225,~
## 11 Prince William Sound (((426007.1 1087250, 426562.5 1088591, 427711.6 1089991~
## 12 Southeast            (((1287777 744574.1, 1290183 745970.8, 1292940 746262.7~
## 13 Yukon                (((-375318 1473998, -373723.9 1473487, -373064.8 147393~
pop <- read.csv("../shapefiles/alaska_population.csv")

head(pop)
##   year     city      lat       lng population
## 1 2015     Adak 51.88000 -176.6581        122
## 2 2015   Akhiok 56.94556 -154.1703         84
## 3 2015 Akiachak 60.90944 -161.4314        562
## 4 2015    Akiak 60.91222 -161.2139        399
## 5 2015   Akutan 54.13556 -165.7731        899
## 6 2015 Alakanuk 62.68889 -164.6153        777
class(pop)
## [1] "data.frame"
pop_4326 <- st_as_sf(pop, 
                  coords = c('lng', 'lat'),
                  crs = 4326, #4236 is WGS84. You have to read in the data in the format it's in and then transform it. You cannot just assign it a diff CRS (such as 3338)
                  remove = F)

head(pop_4326)
## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -176.6581 ymin: 51.88 xmax: -154.1703 ymax: 62.68889
## Geodetic CRS:  WGS 84
##   year     city      lat       lng population                   geometry
## 1 2015     Adak 51.88000 -176.6581        122    POINT (-176.6581 51.88)
## 2 2015   Akhiok 56.94556 -154.1703         84 POINT (-154.1703 56.94556)
## 3 2015 Akiachak 60.90944 -161.4314        562 POINT (-161.4314 60.90944)
## 4 2015    Akiak 60.91222 -161.2139        399 POINT (-161.2139 60.91222)
## 5 2015   Akutan 54.13556 -165.7731        899 POINT (-165.7731 54.13556)
## 6 2015 Alakanuk 62.68889 -164.6153        777 POINT (-164.6153 62.68889)
class(pop_4326)
## [1] "sf"         "data.frame"
#pop_joined <- st_join(pop_4326, ak_regions_3338, join = st_within) comment out bc not working, don't need?
pop_3338 <- st_transform(pop_4326, crs = 3338)

head(pop_3338)
## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -1537925 ymin: 472626.9 xmax: -10340.71 ymax: 1456223
## Projected CRS: NAD83 / Alaska Albers
##   year     city      lat       lng population                   geometry
## 1 2015     Adak 51.88000 -176.6581        122  POINT (-1537925 472626.9)
## 2 2015   Akhiok 56.94556 -154.1703         84 POINT (-10340.71 770998.4)
## 3 2015 Akiachak 60.90944 -161.4314        562  POINT (-400885.5 1236460)
## 4 2015    Akiak 60.91222 -161.2139        399  POINT (-389165.7 1235475)
## 5 2015   Akutan 54.13556 -165.7731        899 POINT (-766425.7 526057.8)
## 6 2015 Alakanuk 62.68889 -164.6153        777  POINT (-539724.9 1456223)
pop_joined <- st_join(pop_3338, ak_regions_3338, join = st_within)

head(pop_joined)
## Simple feature collection with 6 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -1537925 ymin: 472626.9 xmax: -10340.71 ymax: 1456223
## Projected CRS: NAD83 / Alaska Albers
##   year     city      lat       lng population region_id           region
## 1 2015     Adak 51.88000 -176.6581        122         1 Aleutian Islands
## 2 2015   Akhiok 56.94556 -154.1703         84         6           Kodiak
## 3 2015 Akiachak 60.90944 -161.4314        562         8        Kuskokwim
## 4 2015    Akiak 60.91222 -161.2139        399         8        Kuskokwim
## 5 2015   Akutan 54.13556 -165.7731        899         1 Aleutian Islands
## 6 2015 Alakanuk 62.68889 -164.6153        777        13            Yukon
##   mgmt_area                   geometry
## 1         3  POINT (-1537925 472626.9)
## 2         3 POINT (-10340.71 770998.4)
## 3         4  POINT (-400885.5 1236460)
## 4         4  POINT (-389165.7 1235475)
## 5         3 POINT (-766425.7 526057.8)
## 6         4  POINT (-539724.9 1456223)
pop_region <- pop_joined %>% 
  as.data.frame() %>% 
  group_by(region) %>% 
  summarise(total_pop = sum(population))

head(pop_region)
## # A tibble: 6 x 2
##   region           total_pop
##   <chr>                <int>
## 1 Aleutian Islands      8840
## 2 Arctic                8419
## 3 Bristol Bay           6947
## 4 Chignik                311
## 5 Cook Inlet          408254
## 6 Copper River          2294
pop_region_3338 <- left_join(ak_regions_3338, pop_region)
## Joining, by = "region"
#joining by region

plot(pop_region_3338)

pop_mgmt_3338 <- pop_region_3338 %>% 
  group_by(mgmt_area) %>% 
  summarize(total_pop = sum(total_pop), do_union = F) #do union to get rid of holes???

head(pop_mgmt_3338)
## Simple feature collection with 4 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2175328 ymin: 405653 xmax: 1579226 ymax: 2383770
## Projected CRS: NAD83 / Alaska Albers
## # A tibble: 4 x 3
##   mgmt_area total_pop                                                   geometry
##       <dbl>     <int>                                         <MULTIPOLYGON [m]>
## 1         1     67252 (((1287777 744574.1, 1290183 745970.8, 1292940 746262.7, ~
## 2         2    417595 (((561012 1148301, 559393.7 1148169, 557797.7 1148492, 55~
## 3         3     24224 (((-1156666 420855.1, -1159837 417990.3, -1161898 416944.~
## 4         4    137476 (((571289.9 2143072, 569941.5 2142691, 569158.2 2142146, ~
plot(pop_mgmt_3338["mgmt_area"])

ggplot(pop_region_3338) +
  geom_sf(aes(fill = total_pop)) 

ggplot(pop_region_3338) +
  geom_sf(aes(fill = total_pop)) +
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high =  "firebrick", labels = comma)

rivers_3338 <- read_sf("../shapefiles/ak_rivers_simp.shp")
st_crs(rivers_3338)
## Coordinate Reference System:
##   User input: Albers 
##   wkt:
## PROJCRS["Albers",
##     BASEGEOGCRS["GCS_GRS 1980(IUGG, 1980)",
##         DATUM["D_unknown",
##             ELLIPSOID["GRS80",6378137,298.257222101,
##                 LENGTHUNIT["metre",1,
##                     ID["EPSG",9001]]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["Degree",0.0174532925199433]]],
##     CONVERSION["unnamed",
##         METHOD["Albers Equal Area",
##             ID["EPSG",9822]],
##         PARAMETER["Latitude of false origin",50,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-154,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",55,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",65,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]
head(rivers_3338)
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: -622169.8 ymin: 557375.7 xmax: 738652.5 ymax: 2316718
## Projected CRS: Albers
## # A tibble: 6 x 4
##   StrOrder region               n                                       geometry
##      <int> <chr>            <int>                          <MULTILINESTRING [m]>
## 1        3 Aleutian Islands   316 ((-616509.9 557976.9, -616240.8 557375.7), (-~
## 2        3 Arctic           10869 ((-20959.19 2009003, -20804.35 2009054), (-20~
## 3        3 Bristol Bay       4884 ((-231927.5 794449.9, -232004.6 794730.2), (-~
## 4        3 Chignik            141 ((-302874.8 715936.3, -302829.4 715920.8), (-~
## 5        3 Cook Inlet        3447 ((-11441.06 987666.2, -11772.61 988142.6), (-~
## 6        3 Copper River      1784 ((540428.6 1180052, 539410.1 1176802), (54043~
ggplot(pop_region_3338) +
  geom_sf(aes(fill = total_pop)) +
  geom_sf(data = rivers_3338, mapping = aes(size = StrOrder)) +
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high =  "firebrick", labels = comma)

ggplot(pop_region_3338) +
  geom_sf(aes(fill = total_pop)) +
  geom_sf(data = rivers_3338, mapping = aes(size = StrOrder)) +
  geom_sf(data = pop_3338, mapping = aes(), size = 0.5) +
  scale_size(range = c(0.01, 0.2), guide = "none") +
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high =  "firebrick", labels = comma)

pop_3857 <- pop_3338 %>%
  st_transform(crs = 3857)
# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), 
                       c("ymin", "xmin", "ymax", "xmax"))
  
  # Coonvert the bbox to an sf polygon, transform it to 3857, 
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
  
  # Overwrite the bbox of the ggmap object with the transformed coordinates 
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}
bbox <- c(-170, 52, -130, 68)   # This is roughly southern Alaska
ak_map <- get_stamenmap(bbox, zoom = 4)
## Source : http://tile.stamen.com/terrain/4/0/3.png
## Source : http://tile.stamen.com/terrain/4/1/3.png
## Source : http://tile.stamen.com/terrain/4/2/3.png
## Source : http://tile.stamen.com/terrain/4/0/4.png
## Source : http://tile.stamen.com/terrain/4/1/4.png
## Source : http://tile.stamen.com/terrain/4/2/4.png
## Source : http://tile.stamen.com/terrain/4/0/5.png
## Source : http://tile.stamen.com/terrain/4/1/5.png
## Source : http://tile.stamen.com/terrain/4/2/5.png
ak_map_3857 <- ggmap_bbox_to_3857(ak_map)
ggmap(ak_map_3857) + 
  geom_sf(data = pop_3857, aes(color = population), inherit.aes = F) +
  scale_color_continuous(low = "khaki", high =  "firebrick", labels = comma)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Visualize SF objects with leaflet

epsg3338 <- leaflet::leafletCRS(
  crsClass = "L.Proj.CRS",
  code = "EPSG:3338",
  proj4def =  "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
  resolutions = 2^(16:7))
st_crs(pop_region_3338)
## Coordinate Reference System:
##   User input: EPSG:3338 
##   wkt:
## PROJCRS["NAD83 / Alaska Albers",
##     BASEGEOGCRS["NAD83",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["Alaska Albers (meters)",
##         METHOD["Albers Equal Area",
##             ID["EPSG",9822]],
##         PARAMETER["Latitude of false origin",50,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-154,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",55,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",65,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Topographic mapping (small scale)."],
##         AREA["United States (USA) - Alaska."],
##         BBOX[51.3,172.42,71.4,-129.99]],
##     ID["EPSG",3338]]
pop_region_4326 <- pop_region_3338 %>% st_transform(crs = 4326)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = "gray",
                    weight = 1)

m
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)

m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),
                    weight = 1,
                    color = "black",
                    fillOpacity = 1,
                    label = ~region) %>% 
        addLegend(position = "bottomleft",
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")

m
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)

m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),
                    weight = 1,
                    color = "black",
                    fillOpacity = 1) %>% 
        addCircleMarkers(data = pop_4326,
                         lat = ~lat,
                         lng = ~lng,
                         radius = ~log(population/500), # arbitrary scaling
                         fillColor = "gray",
                         fillOpacity = 1,
                         weight = 0.25,
                         color = "black",
                         label = ~paste0(pop_4326$city, ", population ", comma(pop_4326$population))) %>%
        addLegend(position = "bottomleft",
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")

m